Generalized Compressibility in a Glass Forming Liquid 
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We introduce a new quantity to probe the glass transition. This quantity is a linear general- 
ized compressibility which depends solely on the positions of the particles. We have performed a 
molecular dynamics simulation on a glass forming liquid consisting of a two component mixture of 
soft spheres in three dimensions. As the temperature is lowered, the generalized compressibility 
drops sharply at the glass transition. Our results are consistent with the kinetic view of the glass 
transition, but not with an underlying second order phase transition. 
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The glass transition is still not well understood despite 
extensive study. Experimentally the glass transition oc- 
curs when the relaxation time exceeds the measurement 
time and particle motion appears to be arrested. It is 
characterized by both kinetic and thermodynamic fea- 
tures. In the supercooled liquid kinetic quantities such as 
the viscosity and relaxation time grow rapidly as the tem- 
perature is lowered. The glass transition is also reflected 
in thermodynamic quantities, e.g., the specific heat at 
constant pressure has a step-like form and the dielectric 
constant has a peak at a frequency dependent tempera- 
ture. 

Theoretically there have been two main approaches to 
the problem |^,^: dynamic and thermodynamic. The 
first category has been dominated by mode coupling the- 
ory (MCT) in which ideally the relaxation time diverges 
at a temperature Tc above the experimental glass tran- 
sition [|. The kinetic view has produced interesting and 
fruitful concepts such as the influence of the energy land- 
scape on relaxation processes 1^,^ and dynamic inhomo- 
geneities j6|-|^ . The thermodynamic viewpoint attributes 
the glass transition to an underlying phase transition hid- 
den from direct experimental observation by extremely 
long relaxation times ||l|,§,| - 11 1 . In most scenarios there 
is an underlying second order phase transition associated 
with a growing correlation length which produces diverg- 
ing relaxation times as well as diverging static susceptib- 
lities |l|-|l3]. More recently Mezard and Parisi have 
argued that the underlying transition is actually a ran- 
dom first order transition signaled by a discontinuity in 
the specific heat. 

In an effort to better characterize the glass transition 
we introduce a novel probe which we call a generalized 
compressibility . Unlike the specific heat which mon- 
itors energy fluctuations, this linear compressibility is a 
function of the microscopic structure of the system: it 
depends solely on the positions of the particles and not 
on their previous history. It is a thermodynamic quan- 
tity in the sense that it is purely a function of the mi- 
crostate of the system dictated by its location in phase 
space. It is easy to compute numerically, and it is sim- 



pler than the dielectric constant which involves both the 
translation and orientation of electric dipoles. By per- 
forming a molecular dynamics simulation of a two com- 
ponent system of soft spheres, we find that the linear 
generalized compressibility drops sharply as the temper- 
ature decreases below the glass transition temperature 
Tg. The drop becomes more and more abrupt as the 
measurement time increases. 

We now derive expressions for the linear and nonlinear 
generalized compressibility. To probe the density fluc- 
tuations, we follow the approach of linear response the- 
ory and consider applying an external potential ^^(t>{r) 

which couples to the local density p{r) — J2iLi <5(r — fi) 
where denotes the position of the i*'' particle, po is the 
average density. AP has units of pressure and sets the 
magnitude of the perturbation. (/)(r) is a dimensionless 
function of position that must be compatible with the pe- 
riodic boundary conditions imposed on the system, i.e., 
it must be continuous across the boundaries, but is oth- 
erwise arbitrary. This adds to the Hamiltonian H of the 
system a term 
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where we have deflned = jy (Pr(j){r)p{r) = '^^4>{ri). 
p^ is the inner product of and p(r), and we can regard 
it as a projection of the density onto a basis function 
4>{r), i.e., p0 =< >. It weights the density fluctua- 
tions according to their spatial position. The application 
of the external potential will induce an average change 
5p^ in p^: 



Sp4> = {P4,)u - {p<t>)u=o 
where the thermal average {P(p)u is given by 
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The partition function Z — Tie^^'^-^^^^ and (3 is the in- 
verse temperature. For small values of AP, this change 
can be calculated using perturbation theory. Up to third 
order in AP, we find 
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where the cumulant averages are 
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with the thermal average (pp = {p^)u=o- The third or- 
der cumulant, eq.(H), is zero in the liquid phase because 
for every configuration there exists an equivalent config- 
uration with the opposite sign of — (p^) and so we will 
not consider this term any further. We can recast eq. (^) 
as a power series in the perturbation AP: 
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where ks is Boltzmann's constant. In the remainder of 
this paper we will focus our attention on the linear (xi) 
and nonlinear (xni) dimensionless generalized compress- 
ibilities defined by the above expressions. We now discuss 
the choice of the function (j>. We consider applying the 
potential along the direction fi of one of the coordinate 
axes so that (j>{f) = 0(r^). A natural candidate for 0(r^) 
is cos(fc^,r'') with k — 2nn/L, where n = 1,2, ... and 
is the volume V. In this case, p0 is the fc*'' mode of the 
cosine transform of the density. However, it is sufficient 
to consider the simpler function 0(r'^) = jr^'l/L. The 
absolute value means that all the particles feel a force 
along the /ith direction pointing towards the origin. The 
results are very similar to (pir^) = cos{k^r^) for small k 
at a fraction of the computational cost. So our results in 
this paper correspond to = \r^\/ L. This is rather 
like a center of mass. Since the system is isotropic, we 
average over the direction p. 

We have performed a molecular dynamics simulation 
on a glass forming liquid [p^pc| ] consisting of a 50:50 
binary mixture of soft spheres in three dimensions. The 
two types of spheres, labelled A and B, differ only in their 
sizes. The interaction between two particles a distance r 
apart is given by Vapij) = eY^Gap/rf^ + Xaf3{r)] where 
the interaction length aap — {(^a + '^p)/'^, ^b/^a = 1-4 
(a, (3 = A, B). For numerical efficiency, we set the cut- 
off function X^pir) = r/cr^^ - A with A = 13/12i2/i3. 
The interaction is cutoff at the minimum of the poten- 
tial Vaf3{r). Energy and length are measured in units 
of e and a a, respectively. Temperature is given in units 



of e/fcs, and time is in units of OAyiri/e where m, the 
mass of the particles, is set to unity. The equations of 
motion were integrated using the leapfrog method pl| ] 
with a time step of 0.005. During each run the tempera- 
ture was kept constant using a constraint algorithm pl| . 
N = Na + Nb is the total number of particles. The sys- 
tem occupies a cube with dimensions (± L/2, ± L/2, ± 
L/2) and periodic boundary conditions. Since N and L 
are fixed in any given run, the density po — N/ is also 
fixed. 
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FIG. 1. Linear generalized compressibility and diffusion 
constant of a one component system as a function of tem- 
perature for 512 soft spheres. Crystallization is clearly seen 
around T=0.57 (po = 1.1). The measurement time was 10^ 
md steps for each temperature. The compressibility was av- 
eraged over 5 runs while the diffusion is shown for a single 
run. 

As a point of reference we determined the mode cou- 
pling Tc by fitting the data for the relaxation time t{T) 
which is defined as the time when the self part of the in- 
termediate scattering function Fs^a{k,t) falls to 1/e [ po[ . 




FsAk,t) = — 
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where the subscript a refers to the particle type, A or 
B. We choose B particles. fi{t) is the position of particle 
i at time t, and (...) refers to an average over different 
configurations. The wavevector k = 2iTq/L where q is 
a vector of integers. For an isotropic system Fs,a{k,t) 
depends only on the magnitude k — \k\. We choose 
k = kmax = '^T^Qmax/ L where kmax IS the position of the 
first maximum of the partial static structure factor S{k). 
For B type particles we use Qmax — 8.3666. We fit the 
ideal MCT form @ r(T) = A'{T - Tc)'^ with A' = 69, 
7 = 1.6, Tc = 0.306 ± 0.005. In finding Tc, we used 
data from 7 temperatures in the range 0.33 < T < 0.39 
below the caging temperature (Tcago = 0.4) where the 
intermediate scattering function first begins to show a 



2 



plateau. Our fit is consistent with the mode coupling 
theory requirement that 7 > 1.6. 

Our procedure for doing runs is as follows. We start 
each run at a high temperature (T=1.5) and lower the 
temperature in steps of AT = 0.05. At each temperature 
we equilibrate for 10^ molecular dynamics (md) steps and 
then measure the quantities of interest for Nt additional 
steps where = 10^ 2 x 10^ 10^ 3 x 10^ or lO'^. All the 
particles move at each md step. The results are then av- 
eraged over up to 40 different initial conditions (different 
initial positions and velocities of the spheres). 

As a check on our procedure for measuring xi and Xni , 
we consider first the case of the crystal. To this end, we 
consider a system of 512 identical {aA = o'b) particles at 
a density po = 1 • 1 . Figure |^ shows the linear compress- 
ibility and the diffusion constant as a function of temper- 
ature for the single component liquid. At high tempera- 
tures xi has a small slope which becomes steeper at low 
temperature. The salient feature is the very sharp drop 
around T=0.57 of the linear generalized compressibility 
and the diffusion constant. The specific heat, which is 
not shown, has a sharp delta function-like peak around 
T=0.57. The low temperature phase (T<0.57) is a crys- 
tal with sharp Bragg peaks in the structure factor. Upon 
heating and cooling, the transition shows hysteresis. All 
these observations are consistent with the fact that crys- 
tallization is a first order transition. Not shown here 
is the nonlinear compressibility which is zero within our 
numerical error. 

We now examine the response of the two component 
glass forming liquid. All the runs were done at a den- 
sity of Po — 0.6 and as/o'A = 1.4. For these parameters 
crystallization is avoided upon cooling. 
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2. Linear generalized compressibility as a function of 



temperature for different measuring times Im: W (A, 40 
runs), 2 x lO'' (o, 32 runs), 10^ (□, 10 runs), 3 x 10*' (O, 6 
runs) and 10^ (*, 4 runs) md steps. System size is 512 parti- 
cles. Po = 0.6 and ob/c^a ~ 1-4. Inset: T > To subset of the 
same data scaled as described in text. 



Figure g shows the linear generalized compressibility 
as a function of temperature for different run times. The 
compressibility at high temperatures is independent of T 
and about an order of magnitude larger than that of the 
single component fluid. In the vicinity of the glass tran- 
sition xi drops. Notice that as the measuring time <m 
increases, the temperature of the drop decreases and be- 
comes more abrupt. The linear compressibility is propor- 
tional to the width of the distribution of p^. If we regard 
/90 as a generalized center of mass, then the drop in xi 
corresponds to the sudden narrowing of the distribution 
P{P(t,) and the sudden arrest in the fluctuations of p^. 
This behavior can be quantified using a scaling ansatz: 
Xii^MyT) = g(p = tM /t{T)), where the characteristic 
time has the Vogel-Fulcher form t(T) = exp{A/{T~To)). 
The inset of Figure || shows that the data collapse onto a 
single curve with A = 0.75, To = 0.15. (The data could 
not be fitted using t(T) = tmct{T) = A'{T - TcV as 
suggested by simple mode couphng theories j|].) This 
value of To lies below the mode coupling Tc — 0.306, the 
upper bound of the glass transition temperature. Notice 
in Fig. ^ that the drop associated with 10^ md steps oc- 
curs approximately at the mode coupling Tc, indicating 
that our long runs were in equilibrium down to the mode 
coupling transition temperature. This scaling suggests 
that xi becomes a step function for infinite Im and that 
the drop in the compressibility would become a discon- 
tinuity at infinitely long times. This is consistent with 
a sudden arrest of the motion of the particles in the liq- 
uid which is the kinetic view of the glass transition. The 
abrupt drop also appears to be in agreement with Mezard 
and Parisi's proposal that the glass transition is a first 
order phase transition [ pl] |. However, the drop indicates 
that our simulations are falling out of equilibrium and 
therefore we cannot really tell if there is a true thermo- 
dynamic transition. 



As an independent check of the glass transition tem- 
perature, we have calculated the specific heat which is 
shown in Figure ||. Notice that the temperature of the 
peak in the specific heat agrees with the temperature at 
which xi drops. 
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FIG. 3. Specific iieat C'v at constant volume as a func- 
tion of temperature calculated using energy fluctuations. All 
parameters are the same as in Figure 0. 



The behavior seen in Figure ^ is similar to that seen 
in measurements of the real part of the frequency de- 
pendent dielectric function e'{uj) In that case as 
the frequency decreased, the temperature of the peak in 
decreased and the drop in e'iuj) below the peak 
became more abrupt. By extrapolating their data to 
CO = 0, Menon and Nagel argued that e'(w = 0) 
should diverge at the glass transition, signaling a second 
order phase transition. We have looked for evidence of 
this divergence by examining samples of different sizes 
to see if the linear generalized compressibility increased 
systematically with system size. As shown in Figure ^ we 
find no size dependence and no indication of a diverging 
linear generalized compressibility. We also find no size 
dependence for the specific heat and hence, no evidence 
of a diverging specific heat (not shown) . This is corrobo- 
rated by recent MCT calculations of a molten salt which 
find that e'{uj 0) goes to a finite value as the glass 
transition is approached p3]. 



We have found hysteresis at the glass transition by first 
cooling a system of 512 particles to our lowest tempera- 
ture T — 0.1 and then heating in steps of AT = 0.05. As 
before we equilibrate at each temperature for lO** time 
steps and then measure quantities for an additional 10^ 
time steps. Our results are shown in the inset of Figure 
^. Notice the slight hysteresis with the rise in xi upon 
warming being at a slightly higher temperature than the 
drop in xi upon cooling. This hysteresis is consistent 
with the kinetic arrest of motion. 
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FIG. 4. Linear compressibility as a function of temperature 
for different system sizes: 216 (o, 5 runs), 512 (□, 10 runs), 
and 1000 (O, 5 runs) particles. The measuring time was 10® 
md steps in all cases. Inset: Linear generalized compressibil- 
ity as a function of temperature for system of 512 particles 
upon cooling and heating. The measuring time was 10® md 
steps in both cases. The data was averaged over 10 runs. 
Other parameters for the main figure and the inset are the 
same as in Figure |2| 

We now turn to the case of the nonlinear generalized 
compressibility Xni given by eq. (H). We are motivated 
by the case of spin glasses where the nonlinear magnetic 
compressibility diverges at the spin glass transition while 
the linear compressibility only has a cusp [^,^ . There 
have been a few studies of nonlinear response functions in 
real glasses [|l^,^, but these have not found any diver- 
gences. Our results are consistent with this conclusion. 
In particular we find that the nonlinear generalized com- 
pressibility is zero above and below the glass transition 
temperature, though it does show a glitch at the glass 
transition. There is no systematic increase with system 
size, indicating the absence of a divergence. Because Xni 
is sensitive to the tails of the distribution of p^, one must 
be careful to obtain a good ensemble average in the liquid 
above the glass transition temperature. We have done 
this by doing 32 runs, each involving 200,000 time steps, 
with different initial conditions, stringing them together 
as though they were one long run and then taking the 
appropriate averages. This produces a better ensemble 
average of < p| >^ which enters into xm in eq. (^. xm 
also took longer to equilibrate than xi ■ A plot Xni versus 
run time shows that one needs to run at least 10^ time 
steps before Xni appeared to saturate. 

To summarize, we have introduced a new thermody- 
namic quantity which depends solely on the positions of 
the particles and not on their histories. This quantity 
drops abruptly at the glass transition which is compati- 
ble with a kinetic arrest of motion, but not with an un- 
derlying second order phase transition. This generalized 
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compressibility can be measured experimentally. It can 
be directly measured in colloidal experiments which mon- 
itor the positions of the particles . Measurements of 
the width of the distribution of pq^ the spatial Fourier 
transform of the density, would also give the linear gen- 
eralized compressibility. 
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